2. The Duffing Oscillator

In this notebook we will explore the Duffing Oscillator and attempt to recreate the time traces and phase portraits shown on the Duffing Oscillator Wikipedia page

%matplotlib inline
from matplotlib import pyplot as plt

import desolver as de
import desolver.backend as D

D.set_float_fmt('float64')

Using numpy backend

2.1. Specifying the Dynamical System

Now let’s specify the right hand side of our dynamical system. It should be

\[\ddot x + \delta\dot x + \alpha x + \beta x^3 = \gamma\cos(\omega t)\]

But desolver only works with first order differential equations, thus we must cast this into a first order system before we can solve it. Thus we obtain the following system

\[\begin{split}\begin{array}{l} \frac{\mathrm{d}x}{\mathrm{dt}} = v_x \\ \frac{\mathrm{d}v_x}{\mathrm{dt}} = -\delta v_x - \alpha x - \beta x^3 + \gamma\cos(\omega t) \end{array}\end{split}\]

@de.rhs_prettifier(
    equ_repr="[vx, -𝛿*vx - α*x - β*x**3 + γ*cos(ω*t)]",
    md_repr=r"""
$$
\begin{array}{l}
\frac{\mathrm{d}x}{\mathrm{dt}} = v_x \\
\frac{\mathrm{d}v_x}{\mathrm{dt}} = -\delta v_x - \alpha x - \beta x^3 + \gamma\cos(\omega t)
\end{array}
$$
"""
)
def rhs(t, state, delta, alpha, beta, gamma, omega, **kwargs):
    x,vx = state
    return D.stack([
        vx,
        -delta*vx - alpha*x - beta*x**3 + gamma*D.cos(omega*t)
    ])

print(rhs)
display(rhs)

[vx, -𝛿*vx - α*x - β*x**3 + γ*cos(ω*t)]

\[\begin{split}\begin{array}{l} \frac{\mathrm{d}x}{\mathrm{dt}} = v_x \\ \frac{\mathrm{d}v_x}{\mathrm{dt}} = -\delta v_x - \alpha x - \beta x^3 + \gamma\cos(\omega t) \end{array}\end{split}\]

Let’s specify the initial conditions as well

y_init = D.array([1., 0.])

And now we’re ready to integrate!

2.2. The Numerical Integration

We will use the same constants from Wikipedia as our constants where the forcing amplitude increases and all the other parameters stay constants.

#Let's define the fixed constants

constants = dict(
    delta = 0.3,
    omega = 1.2,
    alpha = -1.0,
    beta  = 1.0
)

# The period of the system
T = 2*D.pi / constants['omega']

# Initial and Final integration times
t0 = 0.0
tf = 40 * T

a = de.OdeSystem(rhs, y0=y_init, dense_output=True, t=(t0, tf), dt=0.01, rtol=1e-12, atol=1e-12, constants={**constants})
a.method = "RK87"

a.reset()
a.constants['gamma'] = 0.2
a.integrate()

2.3. Plotting the State and Phase Portrait

# Times to evaluate the system at
eval_times = D.linspace(18.1, 38.2, 1000)*T

from matplotlib import gridspec

fig = plt.figure(figsize=(20, 4))

gs  = gridspec.GridSpec(1, 2, width_ratios=[3, 1])
ax0 = fig.add_subplot(gs[0])
ax1 = fig.add_subplot(gs[1])
ax1.set_aspect(1)

ax0.plot(eval_times/T, a.sol(eval_times)[:, 0])
ax0.set_xlim(18.1, 38.2)
ax0.set_ylim(0, 1.4)
ax0.set_xlabel(r"$t/T$")
ax0.set_ylabel(r"$x$")
ax0.set_title(r"$\gamma={}$".format(a.constants['gamma']))

ax1.plot(a[18*T:].y[:, 0], a[18*T:].y[:, 1])
ax1.scatter(a[18*T:].y[-1:, 0], a[18*T:].y[-1:, 1], c='red', s=25, zorder=10)
ax1.set_xlim(-1.6, 1.6)
ax1.set_ylim(-1.6, 1.6)
ax1.set_xlabel(r"$x$")
ax1.set_ylabel(r"$\dot x$")
ax1.grid(which='major')
plt.tight_layout()

We can see that this plot looks near identical to this plot.

\[\gamma=0.20\]

2.4. Integrating for Different Values of Gamma

Now let’s try to recreate the other plots with the varying gamma values.

gamma_values = [0.28, 0.29, 0.37, 0.5, 0.65]
integration_results = []

for gamma in gamma_values:
    a.reset()
    a.constants['gamma'] = gamma

    if gamma == 0.5:
        a.tf = 80.*T
        eval_times = D.linspace(18.1, 78.2, 3000)*T
        integer_period_multiples = D.arange(19, 78)*T
    else:
        a.tf = 40.*T
        eval_times = D.linspace(18.1, 38.2, 1000)*T
        integer_period_multiples = D.arange(19, 38)*T

    a.integrate()
    integration_results.append(((eval_times, a.sol(eval_times)), a.sol(integer_period_multiples)))

for gamma, ((state_times, states), period_states) in zip(gamma_values, integration_results):
    fig = plt.figure(figsize=(20, 4))

    gs  = gridspec.GridSpec(1, 2, width_ratios=[3, 1])
    ax0 = fig.add_subplot(gs[0])
    ax1 = fig.add_subplot(gs[1])
    ax1.set_aspect(1)

    ax0.plot(state_times/T, states[:, 0], zorder=9)
    ax0.set_xlim(state_times[0]/T, state_times[-1]/T)
    if gamma < 0.37:
        ax0.set_ylim(0, 1.4)
    else:
        ax0.set_ylim(-1.6, 1.6)
        ax0.axhline(0, color='k', linewidth=0.5)
    ax0.set_xlabel(r"$t/T$")
    ax0.set_ylabel(r"$x$")
    ax0.set_title(r"$\gamma={}$".format(gamma))

    ax1.plot(states[:, 0], states[:, 1])
    ax1.scatter(period_states[:, 0], period_states[:, 1], c='red', s=25, zorder=10)
    ax1.set_xlim(-1.6, 1.6)
    ax1.set_ylim(-1.6, 1.6)
    ax1.set_xlabel(r"$x$")
    ax1.set_ylabel(r"$\dot x$")
    ax1.grid(which='major')
    plt.tight_layout()

And we see that these plots replicate those on Wikipedia very nicely except for the chaotic one which is highly sensitive to the numerical integrator and tolerances.

( link to image ) \(\gamma = 0.28\)

( link to image ) \(\gamma = 0.29\)

( link to image ) \(\gamma = 0.37\)

( link to image ) \(\gamma = 0.50\)

( link to image ) \(\gamma = 0.65\)